clear all;
close all;

% [slice,tis,reps,clow,chigh,isCBF,nIters]
f_template_graymatter = 'myProject/template/rc1rT13D_asl_file_on_cpwall_nrtasl_file_gm_nii.nii';
template = cMRI(f_template_graymatter);

TIS = 2:11;
REPS = 5:1:40;
nIters = 100;
nSlice = 14;


qualityCBF = zeros(length(TIS),length(REPS),nIters);
qualityAAT = zeros(length(TIS),length(REPS),nIters);
%
% CBF_mean = zeros(length(TIS),length(REPS),nSlice);
% CBF_std = zeros(length(TIS),length(REPS),nSlice);
%
% AAT_mean = zeros(length(TIS),length(REPS),nSlice);
% AAT_std = zeros(length(TIS),length(REPS),nSlice);

fileName = 'Divasl_CBF_AAT_Gray';
savedName = [fileName,...
    '_REPs',num2str(max(REPS)),...
    '_TIs',num2str(max(TIS)),...
    '_method_cmf'];
dirPath = 'Results/DiverseTIREP_100/';
savedFile = fullfile(dirPath,savedName);
load(savedFile);
dataFull = obj;
% CBF_DATA_FULL = data(1).dCBF.dMRI;
% AAT_DATA_FULL = data(1).dDT.dMRI;


for tis = 1:length(TIS)
    for reps = 1:length(REPS)
        %         for n = 1:nIters
        fileName = 'Divasl_CBF_AAT_Gray';
        savedName = [fileName,...
            '_REPs',num2str(REPS(reps)),...
            '_TIs',num2str(TIS(tis)),...
            '_method_cmf'];
        dirPath = 'Results/DiverseTIREP_100/';
        savedFile = fullfile(dirPath,savedName);
        load(savedFile);
        data = obj;
        
        for n = 1:nIters
            vr = dataFull(n).dCBF.dMRI;
            v0 = data(n).dCBF.dMRI;
            qualityCBF(tis,reps,n) = norm(vr(:)) /(norm(vr(:)-v0(:))+eps);
            
            
            vr = dataFull(n).dDT.dMRI;
            v0 = data(n).dDT.dMRI;
            qualityAAT(tis,reps,n) = norm(vr(:)) /(norm(vr(:)-v0(:))+eps);
            
            qualityCBF(tis,reps,n) = 20*log10(qualityCBF(tis,reps,n));
            qualityAAT(tis,reps,n) = 20*log10(qualityAAT(tis,reps,n));
        end
        
        
        
        disp(['TI = ', num2str(TIS(tis)),' REPEAT = ',num2str(REPS(reps))]);
    end
end
%%
qualityCBF_median = mean(qualityCBF,3);
qualityCBF_std = std(qualityCBF,0,3);
figure(1),
imagesc(REPS,TIS,qualityCBF_median,[0 20]);colormap jet; axis equal;axis xy
M = 5:0.5:40;
hold on;
plot(M,100./M,'-m','Linewidth',1);axis equal;axis([5,40,2,11])
plot(M,200./M,'-m','Linewidth',2);axis equal;axis([5,40,2,11])
plot(M,300./M,'-m','Linewidth',4);axis equal;axis([5,40,2,11])
legend('T\propto 100','T\propto 200','T\propto 300','Location','SouthWest')
xlabel('Sampling Repetitions R');
ylabel('Sampling Diversity M');
colorbar;
hold off;

figure(2),
imagesc(REPS,TIS,qualityCBF_std,[0 3]);colormap jet; axis equal;axis xy; axis tight;
M = 5:0.5:40;
hold on;
plot(M,100./M,'-m','Linewidth',1);axis equal;axis([5,40,2,11])
plot(M,200./M,'-m','Linewidth',2);axis equal;axis([5,40,2,11])
plot(M,300./M,'-m','Linewidth',4);axis equal;axis([5,40,2,11])
legend('T\propto 100','T\propto 200','T\propto 300','Location','SouthWest')
xlabel('Sampling Repetitions R');
ylabel('Sampling Diversity M');
colorbar;
hold off;

%%
qualityAAT_median = mean(qualityAAT,3);
qualityAAT_std = std(qualityAAT,0,3);
figure(1),imagesc(REPS,TIS,qualityAAT_median,[0 20]);colormap jet; axis equal;axis xy
M = 5:0.5:40;
hold on;
plot(M,100./M,'-m','Linewidth',1);axis equal;axis([5,40,2,11])
plot(M,200./M,'-m','Linewidth',2);axis equal;axis([5,40,2,11])
plot(M,300./M,'-m','Linewidth',4);axis equal;axis([5,40,2,11])
legend('T\propto 100','T\propto 200','T\propto 300','Location','SouthWest')
xlabel('Sampling Repetitions R');
ylabel('Sampling Diversity M');
colorbar;
hold off;

figure(2),
imagesc(REPS,TIS,qualityAAT_std,[0 3]);colormap jet; axis equal;axis xy; axis tight;
M = 5:0.5:40;
hold on;
plot(M,100./M,'-m','Linewidth',1);axis equal;axis([5,40,2,11])
plot(M,200./M,'-m','Linewidth',2);axis equal;axis([5,40,2,11])
plot(M,300./M,'-m','Linewidth',4);axis equal;axis([5,40,2,11])
legend('T\propto 100','T\propto 200','T\propto 300','Location','SouthWest')
xlabel('Sampling Repetitions R');
ylabel('Sampling Diversity M');
colorbar;
hold off;

% figure,plot(M-4,150./M-2+1);axis equal;axis([0.5,36.5,0.5,10.5])
%%
% result_CBF_mean = cMRI(CBF_mean);
% result_CBF_std = cMRI(CBF_std);
%
% result_AAT_mean = cMRI(AAT_mean);
% result_AAT_std = cMRI(AAT_std);
%
% result_CBF_mean.showMRI('SNR of CBF, Average');
% result_CBF_std.showMRI('Std of SNR of CBF');
% result_AAT_mean.showMRI('SNR of AAT, Average');
% result_AAT_std.showMRI('Std of SNR of AAT');
